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ABSTRACT 

We study the dynamical evolution of idealised stellar systems by averaging results from 
many iV-body simulations, each having modest numbers of stars. For isolated systems 
with stars of uniform mass, we discuss aspects of evolution up to the point of core collapse: 
relaxation and its iV-dependence, the evolution of the density profile, the development of 
the velocity dispersion and anisotropy, and the rate of stellar escape. 
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1 INTRODUCTION: THE STATISTICAL STUDY OF A-BODY SYS- 
TEMS 



There have been dramatic improvements in the performance of astrophysical A-body 
simulations since the pioneering work of von Hoerner (1960). Several, e.g. regularisation 
of close encounters, have taken place in the domain of software, i.e. developments of the 
programs used. Many of these are documented in Aarseth (1985), and indeed due to him, 
but others are beginning to have an impact on the field of collisional stellar dynamics, 
including tree codes (Barnes & Hut 1986, McMillan & Aarseth 1993) 

Of comparable importance are improvements in hardware, many of which have been 
exploited quickly into collisional stellar dynamics. Among these are special-purpose devices 
such as Grape-2 and its descendants (Ito et al. 1991), and general-purpose machines such 
as the Cray Series supercomputers. The latter have, for example, been used effectively 
by McMillan, Hut & Makino (1990) in their studies of systems with primordial binaries. 
Even so, dedicated workstations are still competitive (Aarseth & Heggie 1992). 

A further group of hardware devices suitable for stellar dynamics are parallel pro- 
cessors, such as the Connection Machine (Hut & Makino 1989). Other machines which 
have been used for problems in stellar dynamics include the Distributed Array Processor, 
which was used in a study of gravothermal behaviour (Heggie 1989), and transputer arrays 
(Sweatman 1993a). Similar machines with faster (vector) processors are now available, 
and have been used in the project described later in this paper. 

The improvement in computational power over the last three decades can be exploited 
in several different ways. One aim is to study larger systems. This is driven by the fact 
that real globular clusters have many more stars than the largest useful computations, by 
at least an order of magnitude. It is not at all clear that this is necessarily the best way 
of using increased computing power, as the discussion of the following section shows. 

Another aim is simply to carry out calculations more quickly. For example, worksta- 
tion calculations on systems with 2500 stars and a substantial proportion of primordial 
binaries may still require as much as 2000 cpu hours (Aarseth & Heggie 1992). 

A third aim is to exploit the power to produce better results, which means better 
statistics. It is well known that little can be done to improve the accuracy of the positions 
and velocities of the stars: errors in these quantities grow exponentially on a time scale 
which is about t cr /8 (Goodman, Heggie & Hut 1993), where t cr is the crossing time. 
Therefore, on the timescales for interesting dynamical evolution (which are of order 10 2 t cr 
for N = 1000) the positions and velocities of the stars are quite wrong. Nevertheless it 
is an article of faith among practitioners of A-body modelling that the statistical results 
are meaningful (e.g. Aarseth & Lecar 1975). Therefore, when we speak of improving the 
results, we mean improving the quality of statistical results. Furthermore, since statistical 
results are the only useful results from these simulations, it is surprising that in many 
cases we are content with a single measurement of the statistical results of importance. Of 
course, it may be that measurements taken at different times on the same system will allow 
a certain amount of statistical sampling, but it is by no means clear how independent these 
are. It is much safer, as stressed, for example, by McMillan, Hut & Makino (1991), to 
compute several simulations with identically distributed initial parameters, differing only 
in the random numbers used to compute each realisation. 
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This approach has been exploited previously by Casertano (1985), who studied as- 
pects such as the distribution of escape energies, though for systems with 8 stars, and by 
McMillan, Caseranto & Hut (1988), who concentrated on an intensive study of the relax- 
ation timescale in iV-body systems with N in the range 16 < N < 1024. Our emphasis 
is on the wider range of dynamical processes, not excluding relaxation, which are of im- 
portance in astrophysical applications, such as the evolution of the spatial distribution of 
stars, anisotropy, escape, binary activity, and so on. 

This paper continues with a summary of the software and hardware issues which are 
relevant to surveys of the kind that we have attempted. We then discuss a number of 
results for isolated systems in which all stars have the same mass: the evolution of the 
space distribution of the stars and the distribution of velocities (including anisotropy), 
and the rate of escape. In this paper we present and discuss data up to the end of 
core collapse. The next paper in this series continues the discussion well into the post- 
collapse regime, and adds some detailed results on the evolution of binaries, and their 
effects on the system. Subsequent papers will deal with systems of stars in which the 
dynamics is more complicated in various ways: systems with a spectrum of masses, tidally 
truncated systems, and systems affected by mass-loss through stellar evolution. Further 
papers will also consider the effects of rotation, which are less easily modeled by Fokker- 
Planck techniques, and mass segregation. 

One purpose of these studies is to test the results of the valuable paper by Chernoff & 
Weinberg (1990), who studied a similar range of problems (up to the point of core collapse) 
with the Fokker-Planck method. More generally we aim to compare our iV-body results 
with those of a variety of simple models and scaling arguments. It might be argued ab initio 
that such an attempt is either bound to fail or to lack meaning, since these simple models 
usually apply in the limit of large N, and it is not clear whether the small models we study 
lie sufficiently far into the asymptotic regime. But in cases where the iV-body results do 
agree with the simple models, we can have some confidence that the simple models are an 
adequate description of the phenomenon in question and its underlying mechanisms. It 
might also be said that we take the approach of a theorist attempting to devise a model 
for a set of observations: starting from familiar physical principles or models (in this case, 
those of stellar dynamics), our aim is to construct the simplest consistent theory which 
accounts for the observations. We allow ourselves to adjust ill-determined parameters to 
optimise the fit, and include no detailed features which are not required by the data. 

2 A PARALLEL STUDY OF iV-BODY STATISTICS 
2.1 Hardware 

Much of the computational work reported below has been carried out on two different 
computers, both installed at the Edinburgh Parallel Computing Centre. One is the Meiko 
transputer array, the other the "Grand Challenge Machine" , also called maxwell, which is 
an array of 64 i860 processors, installed in the autumn of 1990. Though parallel computers 
are often thought of as being very hard to program, it is quite easy to mount a large number 
of independent iV-body problems, as the only messages passed between processors involve 
simple tasks of assigning input parameters and gathering output. The difficulties are not 
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much greater than those involved in mounting an A-body program on any new single- 
processor hardware. 

The transputer array (ECS, for Edinburgh Computing Surface) contains 400 proces- 
sors, each of which has communications and at least 4Mbytes of local storage. As the 
machine is configured for multi-user operation, generally only smaller "domains", of up 
to 131 transputers, are available to each user. Programs are developed and run under 
a software environment known as "CS-Tools". Individual processors are rated at about 
lMflop. In practice, the code NBODY1 (Aarseth 1985) runs at about one fifth times the 
speed obtained on a Sun ELC. 

The Grand Challenge Machine has a similar structure and environment, which made 
the task of porting code to this machine very quick and painless. It differs mainly in the 
speed of the processors, which are roughly 100 times faster than those of the ECS. Each 
i860 is a vector processor, with communications and local storage, which is capable in 
principle of a performance of 10 2 Mflops, and the entire machine has a peak performance 
of about 5Gflops. It is also very cost-effective, at around $1M. The main purpose of the 
machine is computations in high-energy physics, and performance approaching the peak 
has been obtained in this application. But such performance is only available with suitable 
hand-crafted code, and much more modest performance is obtained in routine FORTRAN 
programs; NBODY1 runs at a speed about 5 times faster than on the ECS. On the other 
hand, while it was being installed and tested (in the winter of 1990/91), much computing 
time was available, and was used for some of the calculations described below (see §2.4). 

Later in the course of this project a DEC Alpha superscalar machine (called "fringe" ) 
was installed on field test at the University of Edinburgh. This is a single-processor ma- 
chine, but very fast, and data from large numbers of A-body simulations were built up 
sequentially. 

2.2 Statistical considerations 

As already mentioned in §1, there are at least two ways in which the power of parallel 
computers can be exploited, either in larger simulations or in more numerous simulations. 
In order to assess the value of these two approaches let us consider the i860 machine, with 
64 processors, and let us suppose that it is feasible to study a 1000-body problem on a 
single processor. 

First we consider what can be achieved by using the entire machine for a single large 
calculation. With 64 processors it would be feasible to study a problem in which the 
computational effort was 64 times larger (provided that the system could be used at 100% 
efficiency). With a simple direct summation algorithm, however, the computational effort 
grows roughly as A 3 (Aarseth 1985, with an extra power of N because we are concerned 
with processes occurring on the relaxation time scale), and so it is feasible actually to com- 
pute a 4000-body system. Assuming that statistical results yield estimates with standard 
deviation proportional to A 1 / 2 , this means that the signal-to-noise ratio for these results 
is improved by a factor of 2, compared with what can be achieved in a comparable time 
with a single processor. 

Now suppose the same machine is used to run 64 simultaneous (but independent) 
1000-body simulations (and it is very easy to do so at virtually 100% efficiency). Then 
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the signal-to-noise ratio is improved by a factor of 8 compared with results for a single 
simulation. Therefore it is clear that the statistical quality of the results is optimised by 
adopting the second strategy. The exact figures, and especially the assumed number of 
bodies for a feasible single-processor simulation, are irrelevant. It is easy to see, assuming 
100% efficiency, that the second strategy is better provided that the computational effort 
varies as iV Q for any a > 1. For Ahmad-Cohen schemes the empirical value of a is of order 
2.6 (Aarseth 1985, modified as before), while a tree-based code yields a > 2. 

Before we conclude that it is always better to run large numbers of modest simulations 
than one large one, we should consider other issues than simply statistical quality. There 
are some phenomena which may well behave very differently in small and large systems - 
gravothermal effects, which manifest themselves only for systems containing many thou- 
sands of particle (Goodman 1987) are one. Another phenomenon which would be rather 
difficult to study by the use of large numbers of small simulations is the phenomenon of 
core motions (Makino & Sugimoto 1987, Heggie 1989, Sweatman 1993b), because such 
motions could not be expected to occur coherently in different systems. The statistical 
analysis required would be rather different from the kind that is considered in this paper. 
Possibly the estimation of autocorrelation functions from a large number of simulations 
would be a fruitful approach, but it is not pursued here. Nevertheless, for many phenom- 
ena of interest, such as those studied in this paper, there is no known qualitative difference 
between large and small systems, and the most scientifically productive way of exploiting 
parallelism is to run large numbers of simulations in parallel. 

2.3 Software 

In this project some pilot sets of calculations were computed using a version of 
Aarseth's code NBODY1 (Aarseth 1985), which uses individual time steps. This code 
is quite effective up to the point where binary activity becomes a problem, and was used 
in a preliminary study of core collapse. Subsequently, a code including regularisation of 
binary encounters was adopted, along the lines described briefly in Heggie (1973), but 
extended to include "freezing" of hard isolated binaries (Aarseth 1985). This allowed the 
computation of new series of isolated models, with 250, 500 and 1000 equal-mass stars 
in each, which included the end of core collapse and subsequent reexpansion. This code 
was still too inefficient to deal with very hard interactions between stars and binaries and 
binaries themselves during the advanced collapse and post-collapse phases. Some cases 
were stopped because of the occurrence of excessive energy errors, and this resulted in 
substantial deterioration of overall statistics. Therefore our definitive sets of calculations 
of isolated systems with 250, 500, 1000 and 2000 equal-mass stars were carried out with 
Aarseth's code NBODY5. It was also adopted, with some modifications, for calculations 
including tidal effects, mass loss through stellar evolution, a mass spectrum, etc., which 
will be described in later papers in this series. 

One general point to made about software in a parallel environment is reliability. On 
a scalar machine it is possible, if irksome, to give individual attention to a model when 
some unforeseen combination of circumstances, or a programming error, leads to a serious 
degradation in accuracy. This is not practical if many models are computed simultaneously. 
Then it is necessary to write into the program instructions for automatically restarting a 
computation from a previously stored configuration, but with new parameters which will 
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allow a more accurate integration (McMillan et al. 1991). 

For the organisation of parallel computations some additional programming consider- 
ations are relevant. The domain of processors is conceptually divided into a single 'master' 
and a number of 'slaves', one for each simulation in the set. The computation is initialised 
when the master issues initial data (e.g. N, and a seed for a random number generator) 
to each slave. Each slave runs a version of an iV-body program. When the slave is ready 
to output analyzed data (which in our calculations was done every time unit) it does so 
to the master, which adds each such line of output to a single results file. Each line is 
identified with the time of the output (in iV-body time units) and the identity of the slave. 
The calculations running on different slaves proceed at different rates, and so the lines are 
not in any definite time-order. Nevertheless, after all slaves have finished their simulations, 
the results file can be analyzed, all lines of output which refer to the same time (but from 
different slaves and simulations) being combined statistically. 

Each slave also occasionally produces a file of data which is sufficient to restart the 
run, if this is required because of a degradation in accuracy, or if the run has been inter- 
rupted for any reason, or has to be scheduled over several production periods (in a shared 
environment) . On the other hand this was not regarded as part of the "scientific" output 
of the project, as these data sets were regularly overwritten. 

2.4 Initial conditions and data 

In this section we discuss the initial parameters of the runs, summarise the data which 
was output, and present some statistics on the sets of runs which have been completed and 
are relevant to the results of this paper. 

The initial conditions for these models were drawn from a Plummer model, with equal 
masses, and no primordial binaries (except for any which occur by accident). The systems 
are isolated, with no external tidal effects. Units, which may be referred to below as 
"AT-body units", are standard (Heggie & Mathieu 1986), and lead to the following initial 
values: the crossing time is 2y/2, the virial radius is 1, and the initial rms 3-dimensional 
speed is 1/ y/2. These units are achieved by first shifting to the barycentric frame and then 
rescaling the positions and velocities so that the initial kinetic and potential energies take 
the required values. 

At each time unit, each model outputted a list of data. For these systems with equal 
masses it was not felt necessary to store data on individual particles, and so most of the 
items are statistical. To some extent the list evolved throughout the course of our project, 
as new aspects of interest emerged. The following list is comprehensive, but for the reason 
just stated, not all of this information is available for all sets of calculations. 

(1) the number of the run in the current series; 

(2) the time in A^-body units; 

(3) the total mass of non-escapers (defined in §3.3); 

(4) the virial ratio (i.e. the ratio of kinetic to binding energy of all bound stars, excluding 
"internal" kinetic and potential energies of recognised binaries); 

(5) the central potential (i.e. the potential at the density centre); 

(6) the core radius r c (defined as in Casertano & Hut (1985)). In addition, the core radius 
was also computed from the definition r\ = 3v 2 /4irGp c , where v c and p c are the central 
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three dimensional velocity dispersion and central density, which in turn were computed 
from information on stars within the first Lagrangian radius (cf.(13) and (31) below); 

(7) the mass of the core, computed in two ways: (i) as in Casertano & Hut (1985), and 
(ii) by summing the masses of stars within a sphere of radius r c (calculated by the second 
method outlined in (6) above) centered on the density centre, which in turn is defined 
below in §3.1.1); 

(8) the ratio of the central density (defined as in Casertano & Hut (1985)) to the mean 
density within the half-mass radius; 

(9) the number of stars in the core (computed in two ways, as in (7) above); 

(10) the three-dimensional root mean square speed of stars in the core (except that a 
binary, defined here as a regularised bound pair, contributes only the speed of motion of 
its bary centre); 

(11) the coordinates and velocity components of the density centre, the latter being defined 
by the change in coordinates over one time unit; 

(12) the coordinates and velocities of the centre of mass; 

(13) Lagrangian radii (radii of spheres containing fixed fractions of the mass of bound 
stars, relative to the density centre) for fractions 1, 2, 5, 10, 20, 30, 40, 50, 75 and 90%; 

(14) the number of bound stars; 

(15) the number of escaping single stars; 

(16) the number of regularised binaries (the parameters which determine whether a pair 
is regularised taking standard values (Aarseth 1985); sometimes we loosely refer to such 
binaries as "hard"); 

(17) the number of unregularised binaries with energy greater roughly than 0.1/cT, where 
kT is defined as 2/3 of the current mean kinetic energy of the single stars and of the 
barycentres of binary stars) ; 

(18) the number of merged pairs (defined as in Aarseth 1985 with standard parameters); 

(19) the number of triple or quadruple configurations (defined as in (18) above); 

(20) the total energy of escaping single stars; 

(21) the total and maximum internal binding energy of all unregularised binaries with 
energies above 0.1/cT; 

(22) the total and maximum internal binding energy of all regularised binaries; 

(23) the total energy of the centres of mass of escaping regularised binaries; 

(24) the total internal energy of escaping regularised binaries; 

(25) the internal binding energy of triple and quadruple configurations (defined as in (19) 
above) ; 

(26) the internal binding energy of mergers (defined as in NBODY5); 

(27) the total internal and external energy of the bound members of the system, defined 
as in Aarseth & Heggie (1992); 

(28) the cumulative integration error in the total energy; 

(29) the maximum and minimum internal energy and radius (relative to the density centre) 
of all regularised binaries; 

(30) the names of stars forming the hardest binary; 

(31) mean square radial and tangential velocities of stars between successive Lagrangian 
radii, and the same quantities corrected for escaping stars which have not yet been removed 
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from the iV-body calculation. 

The output from the various models in each series were collected in a single output 
file. For each series the typical size of an output file is typically 35Mbytes. The file for 
each series was later processed so as to produce statistics for all models in that series at 
each time. As each line is read, the time to which it refers is noted, and the data from 
that line are added to the statistical data for that time. 

Table 1 summarises both the extent of the computational effort in these calculations 
and the numbers of particles and simulations. Smaller numbers of simulations with N > 
1000 have also been included for comparative purposes. Further discussion and results, 
including some data from a 10000-body simulation, are considered in a paper by Giersz 
& Spurzem (1993). It is worth pointing out that the statistical value of these various 
sets of simulations are reasonably well balanced: for example, one set consisting of 40 
simulations with 250 bodies is comparable with one simulation containing 10 stars, as the 
total number of particles is the same. 

Before discussing the results of these simulations we mention here some data about 
their accuracy, as measured by energy conservation. Throughout the core collapse period 
the mean of the absolute value of the total energy error is monotonically increasing and 
reaches at the time of collapse values between 10 -5 to 10 -3 8 for different N. As can be 
expected from simple theoretical arguments (e.g. the larger time to core collapse for larger 
N) the smaller value is for N = 250 and the larger for N = 2000. The energy error for 
individual cases never exceeded 0.002. 

3 STATISTICAL RESULTS FROM iV-BODY SIMULATIONS 

3.1 Spatial evolution 

3.1.1 Evolution of Lagrangian radii 

We now discuss the results from these computations up to about the time of core 
collapse, and begin here with the distribution of mass. We assume spherical symmetry 
relative to the "density centre" (computed as in Aarseth & Heggie 1992). In what follows, 
"Lagrangian radii" are the radii of imaginary spheres containing a fixed fraction / of the N 
bound stars in the system, i.e. excluding escaping stars. Where this fraction corresponds 
to an integral number fN of stars, the radius is that of the fNth star in order of distance 
from the density centre. Where fN is fractional, the value is interpolated. 

Combining data from a large number of runs may be done in several ways, and it is 
not clear at first how consistent the different results will be, or what is the most effective 
method. One could, for example, combine the data on individual stars from all cases 
at a given time, and then compute the Lagrangian radii for this enlarged ensemble of 
stars. Because of the limited data which we elected to output, however, this method could 
not be adopted, and we chose instead a suitable statistical estimator based on the set of 
Lagrangian radii obtained for each case individually. For example, Fig.l shows data for the 
Lagrangian radius defined by the innermost 5% of the mass, for N = 250. These results 
indicate that the mean and median (over all cases computed) behave very similarly, and 
from now on we shall concentrate on the mean. 

The improvement in the quality of the data obtained by combining results from many 
cases in this way is illustrated in Fig. 2, which show the means for all Lagrangian radii 
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in one set of models, Fig. 2a, and the results for a single, typical model, Fig. 2b. The 
improvement in "signal-to-noise" is comparable to what one would expect from simple 
"root n" considerations. It is this improvement which will allow us to re-examine several 
quite old problems in collisional stellar dynamics which have lain relatively fallow for want 
of results of sufficient statistical quality. 

3.1.2 A-dependence of dynamical evolution 

The effect of the number of stars on the evolution of the Lagrangian radii is illustrated 
in Fig. 3. The most obvious feature of this result is the increase of the time scale for core 
evolution with A. This is consistent with the assumption that the evolution is driven 
by two-body relaxation, and that the relaxation time increases nearly linearly with A 
(cf. Spitzer 1987). More precisely, the theoretical result is that the relaxation time scale 
varies as N/ln^N), where 7 is a constant. In what follows we analyse this question 
quantitatively in a number of ways, though the order in which these are described differs 
from the order in which they were carried out. 

One procedural point has to be mentioned first. Fig. 4 shows, for A = 250, a com- 
parison between the iV-body results and two approximate (continuum) models, which are 
based on the theory of relaxation and described more precisely below. The point we wish 
to emphasise is the fact that the A-body results lie initially below those of the contin- 
uum models, which correctly give the value for a Plummer model. The reason for this 
discrepancy can be traced to several sources, some of which arise naturally from the way 
in which the initial conditions are scaled. For example, if 250 point masses are spatially 
distributed as in a Plummer model, the resulting potential will not be exactly that of the 
underlying Plummer model, and it is probable that the mean potential energy of a large 
number of simulations is biased away from its true value. Therefore, when the positions 
of the particles are scaled to the desired value of the total potential, there is a resulting 
bias in the spatial distribution. Another bias arises from use of the potential centre, which 
tends to be a point of overdensity in the system, and therefore Lagrangian radii tend to 
be smaller than they should be. The different values of the initial 2% Lagrangian radii for 
various A are clearly visible on Fig. 3. 

Approximation correction for these biases has been carried out by adding a constant 
correction to each mean Lagrangian radius such as to bring the numerical data and the 
theoretical radius into optimal agreement at time t = 0, or for the average of the values at 
t = and 1, or for the first 10% of the collapse phase. (In the latter two cases allowance 
was made for the evolution of the radii during this relatively short time interval.) It 
may be expected, from our discussion of the sources of these effects, that they are A- 
dependent, and indeed we find that the correction to the innermost Lagrangian radius 
generally decreases with increasing A, and is smaller (in absolute terms) at larger radii. In 
what follows we shall use these adjusted Lagrangian radii for a comparison of the various 
models, though we continue to use Fig. 4 to illustrate the discussion. 

Now we return to the question of how the time scale for evolution varies with A. Fig. 5 
illustrates, for two particular values of A, one way in which this can be done. It shows 
one curve for each Lagrangian radius, and was constructed as follows. At each value of t 
in the 500-body models (which is the abscissa in Fig. 5) the mean Lagrangian radius was 
computed, and then by interpolation we determined the corresponding time at which this 
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mean value of the radius was reached in the 2000-body models. (In fact there may be more 
than one such time, because statistical fluctuations result in evolution of the radius which 
is not quite monotonic. We chose the first such time, and checked that our results are 
insensitive to this choice.) According to the theory of relaxation the ratio of these times, 
which is the scale factor plotted in Fig. 5, should be Sf = N 2 In^N^/Ni ln^A^), where 
in this case N± = 500 and N 2 = 2000. 

The results for all combinations of N± and N 2 show remarkable consistency. The 
values of Sf for the innermost Lagrangian radii (1 - 10%) and outermost Lagrangian radii 
(75 and 90%) are very close together. On the other hand the 20 - 40% Lagrangian radii lie 
below and the 50% radius lies above; this means that these intermediate radii evolve too 
fast in models with smaller N. The explanation for this may be connected with the fact 
that binaries start to influence the evolution relatively earlier in core collapse in low-iV 
systems than for larger N (cf. Giersz & Spurzem 1993 and also Paper 2 in this series). 
The effect of binaries would not be confined to intermediate radii, but that is where their 
influence might be most noticeable in this kind of analysis, because these radii evolve less 
quickly than those further in and further out. 

Fig. 5 shows clear evidence of the Coulomb logarithm in the relaxation time: without 
it we would predict Sf = 4 for the relevant values of N± and N 2 . One would like to be 
able to use this data for an empirical determination of 7, but this is very difficult because 
of the weak predicted dependence of Sf on 7: fluctuations in Sf give unacceptably large 
fluctuations in 7. To determine 7, therefore, we adopted a somewhat more indirect method, 
which we now describe. 

Two-body relaxation theory can be used to study the evolution of stellar systems in 
two ways. One is the gas model (Hachisu & Sugimoto 1978, Lynden-Bell & Eggleton 1980) 
and the other, which is usually thought to be more faithful to the picture underlying relax- 
ation theory, is the Fokker-Planck model (cf. Spitzer 1987). Both models, however, share 
some significant simplifying assumptions, including spherical spatial symmetry, isotropy in 
the velocity distribution, and no loss of mass by escape, at least in the simplest formulations 
of these two models. 

In order to use results from these models to evaluate the evolution of an iV-body 
model, it is necessary to relate the time variables in the different models. Conversion 
between Fokker-Planck and iV-body models depends only on the value of the coefficient 7 
in the expression for the relaxation time, i.e. 

0.065(y/3q) 3 

r prnGHn^N) {) 

(cf. Henon 1975, Spitzer 1987), where a is the rms value of each component of velocity, p 
is the stellar mass density, and m is the individual stellar mass. Conversion between gas 
and N-body models requires in addition the specification of a dimensionless conductivity 
coefficient, denoted by C (Lynden-Bell & Eggleton 1980). A number of previous theoretical 
estimates are listed in Table 2. 

If the time variable of a set of iV-body models is denoted by t^, the scaling between 
this variable and the time variable of a continuum model was carried out in a manner 
analogous to the determination of the scale factor 5/ in Fig. 5. For each value of this 
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gives a single estimate for the scalings which bring the models into agreement. A typical 
result is shown in Fig. 6, though since this illustrates scaling with the Fokker-Planck model 
the result has been expressed in terms of the equivalent value of 7. This case illustrates 
a fairly general finding, that the inner radii agree well within themselves, but that the 
half-mass and 75% radii are relatively discrepant. Also noticeable in the inner radii is a 
tendency for the value of 7 to decrease as increases; this is most noticeable in runs with 
the smallest A, and will be discussed further in §3.1.3.. 

We consider first the results for the radii well inside the half-mass radius. As can be 
inferred from Fig. 6, individual estimates of 7 range widely, especially near the beginning of 
the evolution. Elsewhere, however, the range of values is smaller, and roughly independent 
of time and radius. Examination of results for other sets of calculations also show that it 
is roughly independent of A. The great bulk of the values lie in the range 0.07 < 7 < 0.14, 
and show greater consistency for larger values of A than for A = 250. For the gas model 
the scale factor shows similar fluctuations, but their interpretation is more complicated, 
as the scaling depends on two parameters, viz. 7 and C. 

These data can also be used to estimate global values (one for each value of A) for 
the scalings between the A-body models and the two continuum models. The results were 
then used to estimate the A-independent parameters 7 and C, with results shown in Table 
2. For the gas models, the estimate of both parameters can be made from each pair of 
values of A, which was done in order to investigate the variation of these parameters. The 
values presented in Table 2 are, however, obtained from a least squares fit to the results 
for all A. The Fokker-Planck data yield one estimate of 7 for each A, but again the value 
in Table 2 is a mean. (Note, however, that the process is iterative, as the determination of 
the adjustment to the Lagrangian radii already requires a knowledge of the scaling between 
the A-body and Fokker-Planck models, i.e. a preliminary value of 7.) 

The variations about these means are quite modest. For the Fokker-Planck determi- 
nation of 7, results for the four values of A range from about 0.108 to about 0.120. For the 
gas model the range of values of C is also quite small, from about 0.099 to about 0.106, but 
the range for 7 is greater: 0.104 to 0.138. Note that these variations are correlated, larger 
values of 7 corresponding to (and partially compensating) smaller values of C. Despite 
these variations, the consistency is unexpected, especially when it is considered that 7 is a 
coefficient in the argument of a logarithm which is usually regarded as "slowly varying" . 
An indication of the overall success of the resulting parameters, at least for the inner radii, 
is provided by Fig. 7, which for one value of A shows the comparison between the gas, 
Fokker-Planck and A-body models. 

It should be stressed again that our empirically determined values of 7 and C are 
based mainly on data for the innermost Lagrangian radii. 

Also shown in Table 2 are some earlier theoretical estimates of the parameters 7 and 
C. It is interesting to point out that one of Henon's values of 7 stemmed from a close 
reexamination of Chandrasekhar's theory of relaxation, which he undertook in order to 
improve a previous comparison between the results of A-body and Fokker-Planck models 
(Aarseth, Henon & Wielen 1974). It is only slightly greater than our value. 

3.1.3 Differences between A-body and continuum results 

We now discuss the discrepant results in Fig. 6 (and related data), beginning with 
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the half-mass radius. It is well known that the half-mass radius is nearly constant in the 
collapse of a Plummer model, which can be understood in terms of energy conservation and 
the empirical fact that the half-mass radius, expressed as a fraction of the virial radius, is 
not sensitive to the model under consideration (Spitzer 1987). The implication of this near 
constancy, however, is that the scale factor from which 7 is computed is very sensitive to 
fluctuations in the iV-body data. Nevertheless the difference between the behaviour of the 
half-mass radius and that of the inner radii is systematic, in the sense that the rn grows 
too quickly in the continuum models (Fig. 7), a fact which was already noted in a previous 
small-scale statistical study of 100-body systems (Heggie 1991). As was mentioned in 
§3.1.2 similar behaviour is exhibited by the 20 and 75% Lagrangian radii, and this already 
becomes visible soon after the start of the evolution. 

A possible reason for this is the development of anisotropy in the A-body models, 
which is discussed further in §3.2 below and in much detail by Giersz & Spurzem (1993). 
Their finding is that the evolution of Th is not even very well described by any reasonable 
anisotropic gaseous model. For the half mass radius the best agreement with the N-body 
data is provided by the anisotropic Fokker-Planck models of Stodolkiewicz (1982) who 
uses a Monte-Carlo technique. Since the only relevant published data from these models 
is diagrammatic, a quantitative measure of the agreement is not yet possible. 

Another factor which is neglected in the continuum models is the change in mass (M) 
and energy (E) of the system caused by escape (cf. §3.3 below). We have corrected the 
values of r h in the gas models by using the actual values of mass and energy in the A-body 
models, assuming that ru oc M 2 /E. The correction has a noticeable effect, but does not 
substantially improve the agreement between the two models. 

It is easily seen in Fig. 7 that the discrepancy is worse in the gas model than in the 
Fokker-Planck model, the former evolving more slowly beyond the 10% radius. This has 
nothing to do with the A-body models, of course, but is a deficiency in the gas model. In 
effect it means that the value of C must increase slightly with radius, reinforcing the fact 
that the values in Table 2 refer to the innermost radii. 

The other noticeable discrepant feature in Fig. 6 is the slight decrease in the scaling 
factors with time. It is related to a systematic discrepancy between the evolution of the 
continuum and A-body radii which is visible particularly clearly in Fig. 5. This trend is 
qualitatively consistent with two possible explanations, both of which are of dynamical in- 
terest. One is the assumption, discussed by Spitzer (1987), that, for the core, the argument 
of the Coulomb logarithm in eq.(l) should be proportional to A c , the number of stars in 
the core, rather than A (cf. Table 2). We discuss this possibility further below. The other 
possible explanation is the growing activity of binary stars, which, in such relatively small 
models, begins to influence the collapse of the core well before it brings the collapse to a 
close. (This process is not modeled at all in the Fokker-Planck and gas models used for 
this comparison.) We shall discuss this in the next paper in this series, which analyses the 
behaviour of binaries throughout the evolution, including the phase beyond core collapse. 
As a precaution, however, in the determination of C and 7 we have discarded data at times 
when binaries have absorbed more than 1% of the total energy of the system. 

In order to investigate the effect of a variable Coulomb logarithm we have repeated 
the computation of gas models but have chosen for the Coulomb logarithm the expression 
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ln(7-/V c ) if N(r) < N c and ln(7iV(r)) otherwise, where N c is the number of stars in the 
core (Spitzer 1987, p. 149) and N(r) is the number of stars within radius r. The value of 
7 chosen was 0.1 (cf. Table 2), and the computation carried out for N = 250, for which 
the problems with fitting the Lagrangian radii by results of gas computations seemed most 
difficult. The rationale for the choice of Coulomb logarithm starts from the fact that the 
argument of the logarithm should be the ratio of the largest effective impact parameter to 
the 90° deflection distance. The latter, which depends on the velocity dispersion, varies 
rather slowly throughout the system. The former may be taken as being comparable with 
the radius of that part of the system where the density is not much smaller than the local 
value. Within the core this distance will be of order r c , the core radius, and outside the 
core it is of order the radius r itself. In a nearly isothermal halo (outside the core) r oc N(r) 
approximately. 

The result of this calculation showed an overall slowing of the evolution, since the 
modification of the Coulomb logarithm from ln(7iV) lengthens the relaxation time. But 
the effect decreases as one moves to larger radii. Compared with the iV-body models with 
N = 250, it was found that the gas model with variable Coulomb logarithm gave results 
for the scale factor (between time variables for the gas and iV-body models) which were 
a little more consistent for different radii (but still excluding the half-mass radius). On 
the other hand the time-dependence of the scale factor was virtually unimproved, which 
suggests that this has some other explanation, such as binary evolution. 

In conclusion, there is some evidence favoring use of a variable Coulomb logarithm. In 
addition, the comparison between iV-body and Fokker-Planck data gives some indication 
that 7 may increase slightly with increasing radius, which is qualitatively consistent with 
our prescription for a variable Coulomb logarithm . Unfortunately this result cannot be 
conclusive because of the large fluctuations in the values of 7 (cf. Fig. 6). 

3.2 The velocity distribution 

Fig. 8 shows some typical results for the velocity dispersions (radial and tangential, 
or the sum of these). The general increase with time in the inner parts of the cluster 
(Fig. 8a) is to be expected from core collapse, and the results from smaller N models show 
a very similar trend with greater noise. The illustrated comparison with the results of 
gas calculations is satisfactory until late in the core collapse phase, and there is even 
a suggestion that the iV-body models exhibit the initial cooling, which is also found in 
Fokker-Planck models (cf. Cohn 1980, Fig. 7). There is no discernible anisotropy in the 
innermost shell, and slight evidence (seen in a graph of the quantity (v^)/{vf), which 
is not, however, reproduced here) between the Lagrangian radii for 20 and 50% of the 
mass. Anisotropy clearly grows, however, in the outer shells (Fig. 8b). There is a slight 
decline in the radial velocity dispersion, but the tangential dispersion decreases much more 
rapidly. The overall downward trend of the mean velocity dispersion agrees well with the 
predictions of an isotropic gas model and is associated with the general expansion of this 
zone, but the radial dispersion is maintained by the ejection of stars from the inner parts 
of each system. The growth of anisotropy (measured by A = 2 — (vf)/(v^)) with time is 
quite linear, and it reaches values between 0.9 and 1.1 (the larger value corresponding to 
the smallest value of N) at the end of core collapse. Detailed comparisons with anisotropic 
gaseous models are presented in Giersz & Spurzem (1993). 
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3.3 Escape 

3.3.1 Review of previous work 

The rate at which stars escape from an isolated clusters of equal stars is a problem with 
a long history (Table 3). Early theoretical estimates (Ambartsumian 1938, Chandrasekhar 
1943) were based on the theory of relaxation, and led to predictions of a fractional escape 
rate per relaxation time which is independent of N. This translates to an escape rate per 
crossing time which is proportional to log N. Relaxation being understood as a diffusive 
process, Henon (1961) (who also lists some other theoretical estimates of the escape rate) 
then pointed out the difficulty of diffusing stars across the energy of escape, when their 
periods should become extremely long. He therefore gave an estimate of escape from two- 
body interactions, though the result depends on the assumed distribution of velocities; it 
is 

dN 256 v / 2VG 2 m 2 f 2j f f (e + e' - </>) 3 / 2 , 

W = § 1 rdr j j 7 /(£)/(£ )d£de ' (2) 

where m is the individual stellar mass, r is the distance from the centre of the cluster, 
f(e) is the phase-space density expressed as a function of the specific energy e, 0(r) is the 
potential, and the integration is taken over the range e < 0, e' < 0, e + e' > 0. Henon's 
point of view was further modified by Spitzer & Shapiro (1972), who pointed out that the 
distribution itself evolves on a relaxation time scale, and then a star which has diffused to 
energies a little below the escape limit can escape in a single two-body encounter in the 
core. Some previous numerical estimates for the escape rate are also summarised in Table 
3. 

3.3.2 Results from the present iV-body simulations 

Escaping stars in our iV-body models are identified according to the following defini- 
tion. Their energy (computed in the rest frame of the centre of mass of the entire system) 
must be positive, and their distance from the density centre must exceed a boundary radius 
which is chosen as 10 or 20 times the half-mass radius. The results presented in this paper 
are for boundary radius equal to 20 times r^, except where stated in §3.3.3. 

Two typical sets of results are shown (along with a theoretical comparison, which we 
discuss below) in Fig. 9. Note that these are means over many runs. We have checked that 
the mean agrees quite well with the median number of escapers, but there is a wide range 
about the mean. For N = 250, for example (not shown), the mean number of escapers at 
the end of core collapse is approximately 4.0, but in individual models may range from 
to 9. Note the clear increase in the rate of escape with time; such an increase was already 
noted in the Monte Carlo models of Spitzer & Shull (1972), who found that it increased 
by a factor of about four during the course of core collapse. We found that it increased by 
a factor of about 6 for N = 250. For higher N the increase is even larger (Fig. 10). 

For a first comparison between our own results and those of theory we have adopted 
a synthesis of the general theoretical result of Henon with the ideas of Spitzer & Shapiro. 
What we have done is to apply Henon's result to the time-dependent distribution function 
given by the (isotropic) Fokker-Planck model. Results of such a comparison are shown in 
Fig. 9, which shows that the agreement is rather poor, except initially. Up to the end of 
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core collapse the mean number of escapers in the A-body models exceeds the theoretical 
prediction by a factor of about 1.8 and 2.8 for A = 500 and 1000, respectively. The 
results for other values of A are qualitatively similar, and the agreement at early times is 
clearer in the larger models than smaller one (as we can see in Fig. 9b). For all models (but 
particularly for this with low A) the comparison is complicated by the fact that stars are 
only deemed to escape after they have reached a sufficiently large radius. This depresses 
the very early times (cf. also Fig. 10), and makes the number of escapers 

depend on the definition of the boundary radius. 

The disagreement between iV-body and Fokker-Planck results is not surprising, be- 
cause eq.(2) is based on the assumption that the distribution of velocities is isotropic. This 
is true of all our theoretical comparisons, but it is especially critical in this case, because 
escape tends to take place for stars which are already loosely bound and yet have suffi- 
ciently small angular momenta to pass through the core. In addition to anisotropy, which 
we discuss further below, further possible explanations of this discrepancy are explored in 
Paper II. 

Fig. 11 shows the situation for the energy carried off by escapers, i.e. their asymptotic 
kinetic energy at infinity. What is plotted here is, as usual, the mean value over all the 
simulations at this value of A. Use of the median would be particularly unsuitable for 
escaper statistics, as a substantial fraction of cases may have no escapers for a significant 
part of the core collapse phase. 

In Fig. 11 the theoretical result is based on a general expression given by Henon (1969) 
for isotropic distribution functions, and evaluated by him for a Plummer model with un- 
equal masses. As before we applied it to the evolving distribution function produced by 
an isotropic Fokker-Planck code. Again there is agreement at early times, but a faster 
rise throughout much of core collapse. Even more dramatic is the abrupt rise towards the 
close of core collapse, and it is natural to interpret this in terms of escapers produced in 
three-body encounters. Again discussion of this process is relegated to Paper II in this 
series. 

3.3.3 The effect of anisotropy 

In order to understand the influence of anisotropy on the rate of escape it would 
be desirable to calculate the evolution of the distribution function using an anisotropic 
model, and then repeat our computation of the escape rate using a suitable anisotropic 
generalisation of Henon's formula, eq. (2). In fact anisotropic gas models have been 
developed by several authors (Larson 1970, Bettwieser 1983, Louis and Spurzem 1991), 
and a recent version of such a code is compared with the results of our and other A-body 
data in Giersz & Spurzem (1993). Unfortunately, however, the gas model does not directly 
yield the distribution function in phase space, and so these results cannot be applied to 
the discussion of the escape rate. 

Since the effect of anisotropy on escape appears to be quite unknown, we have formu- 
lated a model problem which should provide a semi-quantitative guide. What we have done 
is to use a Monte Carlo method to compute the escape rate in a sequence of anisotropic 
models. The sequence we chose was devised by Dejonghe (1987), and has the benefit that 
all the models in the sequence have the same spatial density distribution as in Plummer's 
model. This allows us to compare our results with Henon's well known results for this case 
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(Henon 1961, 1969). More importantly, the models have the same density distribution as 
our initial conditions; thus the study of this sequence of models allows us to determine 
the effects of anisotropy, while our study of the Fokker-Planck model allows us to study 
the effect of dynamical evolution. Dejonghe's models have one scale-free free parameter, q, 
which we took in the range from (Plummer's model) to 1.5. It is related to the anisotropy 
A defined in §3.2 by 

q = A(l + l/r 2 ). (3) 

In order to compute the escape rate the following procedure was used. First, the radius 
r is selected from a probability density function proportional to 47rr 2 n(r) 2 , where n is the 
number-density. Next, the velocities of two stars are selected from the local distribution of 
velocities. Then the impact parameter p is chosen with probability density proportional to 
p up to a maximum value Pmax, and an angle which determines the relative orientation of 
the two stars at the moment when they would be at their minimum separation during an 
encounter, if their paths were undeflected by it. Then it is possible to determine whether 
either star has a speed, after the encounter, above the local escape speed. Then the escape 
rate is the average value of v, the relative speed of the two stars before the encounter, 
suitably normalised. The resulting value is subject to statistical uncertainty, and depends 
on the choice of p m ax- Experimentation with different values, and different numbers of 
trials, allowed us to determine satisfactory choices. 

The results of our calculations are shown in Figs. 12 and 13. Units are standard, i.e. 
those in which M = 1, G = 1 and E = —1/4. For q = the results may be compared 
with those given by Henon, i.e. A = 0.00942 and E = 0.000741. Evidently the presence of 
anisotropy can lead to a large increase in the rate of escape and the flux of energy carried 
off by escaping stars. 

In applying these results to a comparison with our A-body data we have made two 
rather arbitrary assumptions. One is that the effects of evolution (discussed in the previous 
section) and anisotropy can be combined multiplicatively. The second is that we can 
estimate the effect of anisotropy by using the analytical (Plummer-like) model which has 
a value of q given by eq.(3), in which A is taken to be the average anisotropy in the 
outer half of the cluster and r is the 75% Lagrangian radius. When this is done it is 
found that the predicted number and energy of escapers agree with those from the A- 
body simulations quite well as can be seen in Figs. 9 and 11. The agreement is particularly 
good for A > 1000. For smaller A the predicted number and energy of escapers are larger 
then those obtained from the A-body data, and for A = 250 the discrepancy is about 
40% of the iV-body value. The largest discrepancy that we have noticed is the energy 
flux at late times in core collapse; as already mentioned it rises rather suddenly as the 
end of core collapse is approached. It is probable that the increased contribution from 
three-body escape events (i.e. those associated with binaries; cf. Paper II) is beginning 
to make a significant contribution. Another explanation which we cannot rule out is the 
inappropriateness of one or other of the assumptions mentioned earlier in this paragraph, 
though the abruptness of the rise suggests otherwise. 

As we mentioned before the A-body data depend on the definition of the boundary 
radius. If the boundary radius is reduced to lOr^ the predicted values agree with those 
from A-body models to better than 25% throughout most of core collapse. But now, for 
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all values of N that we have studied, they are smaller than the iV-body values. Therefore 
we can expect that in the case of instantaneous removal of escapers from iV-body models 
(which is in effect what happens in the theoretical model) the discrepancies between com- 
puted and predicted values will be larger than for the model with boundary radius 10r/j,. 
But at any rate, our conclusion is that anisotropy has a large enough effect on the escape 
rate that it may be quite enough to reconcile theoretical predictions with our numerical 
results. 

4 DISCUSSION AND CONCLUSIONS 

We have computed large numbers (~ 200) of isolated iV-body simulations with N in 
the range from 250 to 2000. All stars have equal masses, and initially the systems are in 
dynamic equilibrium. The purpose of these studies is to improve the quantitative results 
of such computations by statistically combining data from many simulations. 

We find that the evolution of the spatial distribution of stars, up to the end of core 
collapse, is quantitatively quite consistent with the theory of relaxation. Indeed we have 
been able to estimate a rather reliable value for the Coulomb logarithm. Small deviations 
from the predictions of isotropic models based on the theory of relaxation are present 
at small radii (where they are probably associated with incipient binary activity) and at 
the 50 and 75% Lagrangian radii (where anisotropy of the velocity dispersion may be 
responsible) . 

Inside the radius containing half the mass of the cluster, the evolution of the velocity 
dispersion is also consistent with simplified models in which the distribution of velocities 
is assumed to be isotropic, but there is strong growth of anisotropy in the outer half of the 
mass. A quantitative comparison with the anisotropy predicted by a simplified theory of 
relaxation is presented elsewhere (Giersz & Spurzem 1993). 

The rate of escape of stars is roughly consistent with a hybrid model based on existing 
theory for escape in evolving, isotropic systems, with a somewhat ad hoc but important 
modification due to anisotropy. The escape is assumed to take place by two-body interac- 
tions from stars whose distribution function slowly changes in accordance with the usual 
theory of relaxation. 
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Table 1 

OUTLINE OF THE COMPUTATIONS 



N number of cases cpu (hrs) machine 

250 56 50 super 

500 56 200 super 

1000 40 190 maxwell 

1000 50 1200 fringe 

2000 16 1350 super 

2500 2 13000 copernicus (Sun ELC) 



Note: For parallel machines (super and maxwell) the total number of processor-hours is 
the product of columns 2 and 3. 
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Table 2 



ESTIMATES OF RELAXATION PARAMETERS 



Parameter Source Notes 

C 7 

0.25 Ambartsumian (1938) 

0.35 Chandrasekhar (1942) 

0.15 Henon (1975) Depends slightly on 

velocity distribution 

0.2 McMillan (pers. comm.) 

0.4 Spitzer (1987) 

2N C /N Spitzer (1987) For the core; N c = 

number of stars in core 

0.104 Heggie & Stevenson (1988), 

Heggie (1985), Goodman (1987) 
0.104 0.11 This project (gas model) 

0.115 This project (F-P model) 
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Table 3 



THE ESCAPE RATE 



ooui ce 


+ rlNlrli 
— Z cr ttl\ 1 (IT 


Ambartsumian(1938) 


0.16 In A 


Chandrasekhar(1943) 


0.33 In A 


Henon(1961) 


0.027 


Spitzer & Shull (1972) 


0.06 In A 


Wielen (1974) 


~ 0.1 


Aarseth & Lecar (1975) 


0.2-0.3 


Stodolkiewicz (1982) 


0.3 


Heggie (1991) 


0.09 



Notes 

Relaxation; In A = Coulomb logarithm 
Relaxation and dynamical friction 
Plummer model; 2-body encounters 
Mean value in core collapse 

50 < N < 250 
Plummer model (N = 250) 
Maen value in core collapse 
N = 100 
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FIGURE CAPTIONS 



Fig.l Statistics of the 5% Lagrangian radius for 56 cases with A = 250. The curves 
correspond to the minimum, median, mean and maximum at each time. 
Fig.2 Mean Lagrangian radii for (a) 56 cases with A = 250, compared with (b) a single 
typical case. 

Fig. 3 Mean Lagrangian radii for all series, for a mass fraction of 2%. For a Plummer 
model the value is 0.166. 

Fig.4 Comparison between A-body, gas (IGM) and Fokker-Planck (F-P) models for the 1% 
Lagrangian radius and A = 250. The time scales of the continuum models have been scaled 
empirically to provide a reasonable global fit; the values used were C = 0.104, 7 = 0.11 
(cf. Table 2). 

Fig. 5 The scale factor in time which brings A = 500 and A = 2000 models into agreement. 
All ten Lagrangian radii are shown. The discrepant ones are 20% — 40% radii (below the 
others) and half-mass radius (above the others). The thick line shows scale factor for 
7 = 0.11. The abscissa is the time variable in the 500-body calculations. 
Fig. 6 The value of 7 which brings Fokker-Planck and iV-body results into agreement at 
each time, for A = 500. The discrepant Lagrangian radii (below the others) are the 50% 
and 75% Lagrangian radii. 

Fig.7 Comparison between the gas (IGM), Fokker-Planck (F-P) and iV-body models with 
7 Lagrangian radii for A = 1000, to indicate the overall success of our fitted values of the 
relaxation parameters 7 and C (cf. Table 2). 

Fig. 8 Velocity dispersions for iV-body models with iV = 1000. (a) Total velocity dispersion 
between Lagrangian radii for 1% and 2% of the mass, compared with a gas model, (b) 
Radial and transverse dispersions between Lagrangian radii for 75% and 90%. 
Fig. 9 Number of escapers as a function of time for different sets of models. F-P-I - 
isotropic Fokker-Planck model, P-A-C - anisotropic Plummer model with evolution (see 
text for discussion), (a) for A = 500, (b) for A = 1000. 

Fig. 10 The escape rate as a function of time for A = 1000. The squares represent data 
given by Spitzer & Shull (1972) for a Monte-Carlo model. Their time unit has been scaled 
to A-body time for A = 1000. The triangle represents data obtained by Stodolkiewicz 
(1982) for a Monte-Carlo model. The dot represents the result given by Henon (1961) for 
a Plummer model. The labels 10r/i and 20r h correspond to the boundary radii chosen in 
different sets of A-body simulations (see text for discussion). 
Fig. 11 Energy of escapers as in Fig. 9. 

Fig. 12 Escape rate for an analytic sequence of anisotropic models (the anisotropy being 
determined by q). Error bars are estimated 1-cr confidence limits. The dashed line and the 
solid curve are linear and quadratic fits of the form a + bq and a + bq + cq 2 respectively. 
The linear fit is adequate for all except the lowest value of q. 

Fig. 13 Rate at which energy is carried off by escapers in a sequence of analytical models. 
For other details see the caption to Fig. 12. 
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